%% VK_KERNEL_COMPUTE Compute a viability kernel approximation
%
% SYNOPSIS
%   This function takes a constraint set, an array of differential
%   equations, and a maximum absolute control magnitude, and attempts to
%   compute an approximate viability kernel, by dividing the state-space
%   into a discretised set of points (according to the `discretisation'
%   option specified in `options'), and calling a viability-determination
%   algorithm (usually vk_viable) against each point.
%
%   In addition to the three arguments that this function accepts, additional
%   options can be passed in either as (name, value) pairs, or as a structure
%   generated by vk_options.
%
%   vk_kernel_compute makes use of cellfun by splitting the problem space into
%   discretisation-many sub-problems, which are then passed into cellfun.
%   This is useful because in GNU Octave, parcellfun can be used as a
%   drop-in replacement for cellfun to simultaneously consider the
%   viability of multiple points, using parallel processing.  For MATLAB(R), we
%   have written an implementation of cellfun called vk_cellfun_parfor which
%   facilitates parallel processing when the Parallel Toolkit is available.
%   The choice of cell function to use can be altered by changing the `cell_fn'
%   option (see vk_options).
%
% USAGE
%   % Standard way of calling:
%   V = vk_kernel_compute(k, f, c)
%
%       - `K' is the constraint set, a row vector twice as long as the number
%         of variables,
%
%   % Passing in an options structure, constructed by vk_options:
%   V = vk_kernel_compute(K, f, c, options)
%
%   % Using the default options, except for some specified here:
%   V = vk_compute(K, f, c, ...
%       'name1', value1, ...
%       'name2', value2 [, ...])
%
%   % Using an options structure, and modifying some parameters:
%   V = vk_kernel_compute(K, f, c, options, ...
%       'name1', value1, ...
%       'name2', value2 [, ...])
%
% EXAMPLES
%   % Compute a simple viability kernel
%   K = [0, 1, 0, 1]    % Two dimensions, each with the same upper and
%                       % lower bounds.
%   f = @(x, u) [1/2*x(1) + x(2)*u; u];
%   V = vk_kernel_compute(K, f, 0.001);
%
%   % Compute the same kernel with a higher discretisation
%   V = vk_kernel_compute(K, f, 0.001, 'discretisation', [50, 50]);
%
%   % Compute the same kernel again, but this time using PARCELLFUN
%   V = vk_kernel_compute(K, f, 0.001, ...
%       'discretisation, [50, 50], ...
%       'cell_fn', @(varargin) parcellfun(2, varargin{:}, 'UniformOutput', false));
%
% Requires:  vk_kernel_compute_recursive, vk_options
%
% See also: cellfun, parcellfun, vk_cellfun_parfor, vk_viable

%%
%  Copyright 2014 Jacek B. Krawczyk and Alastair Pharo
%
%  Licensed under the Apache License, Version 2.0 (the "License");
%  you may not use this file except in compliance with the License.
%  You may obtain a copy of the License at
%
%      http://www.apache.org/licenses/LICENSE-2.0
%
%  Unless required by applicable law or agreed to in writing, software
%  distributed under the License is distributed on an "AS IS" BASIS,
%  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
%  See the License for the specific language governing permissions and
%  limitations under the License.
function [V, NV, viable_paths nonviable_paths] = vk_kernel_compute(K, f, c, varargin)

    %% Build options.
    options = vk_options(K, f, c, varargin{:});

    %% These are the options that we are interested in.
    numvars = options.numvars;
    discretisation = options.discretisation;

    %% Create grids representing the search space
    %   LINSPACE is used to accomplish this.  Discretisation in each dimension
    %   is potentially different.
    axes = cell(numvars, 1);
    for i = 1:numvars
        axes{i} = linspace(K(2*i - 1), ...
            K(2*i), discretisation(i));
    end

    %% Make a cumulative product along the dims
    cp = cumprod([1; discretisation(1:end-1)]);

    %% Compute the total number of points
    num_pts = prod(discretisation);

    %% Create a function for use with CELLFUN
    %   This function wraps the VK_KERNEL_COMPUTE_RECURSIVE helper function
    %   prepopulating all variables except for i and x.  See
    %   VK_KERNEL_COMPUTE_CELLFN, below.
    fn = @(i) vk_kernel_compute_cellfn(i, axes, cp, discretisation, K, f, c, options);

    %% Call ARRAYFN or PARARRAYFUN
    [pts, viable_cells, viable_path_cells] = options.array_fn(fn, 0:num_pts-1);

    %% Build the viability kernel from the results of the CELLFUN call
    viable = cell2mat(viable_cells);

    viable_idx = find(viable);
    nonviable_idx = find(not(viable));

    % Reshaping is necessary for empty matrices
    V = reshape(cell2mat(pts(viable_idx))', [numel(viable_idx), numvars]);
    NV = reshape(cell2mat(pts(nonviable_idx))', [numel(nonviable_idx), numvars]);

    viable_paths = viable_path_cells(viable_idx)';
    nonviable_paths = viable_path_cells(nonviable_idx)';
end

% This function is called on each point under consideration for viability.
function [x, viable, viable_path] = ...
      vk_kernel_compute_cellfn(i, axes, cp, discretisation, K, f, c, options)

  if (options.report_progress)
    options.progress_fn(i);
  end

  if (options.cancel_test && options.cancel_test_fn())
    return;
  end

  idx = mod(floor(i ./ cp), discretisation);
  x = zeros(length(idx), 1);
  for n = 1:length(axes)
    axis = axes{n};
    x(n) = axis(idx(n) + 1);
  end

  if (options.debug)
    disp(i);
    disp(x');
  end

  try
    [viable, viable_path] = options.viable_fn(x, K, f, c, options);

  catch
    exception = lasterror();
    warning(['Error computing viability of point ', ...
             mat2str(transpose(x)), ...
             ': ', ...
             exception.message]);
    viable = false;
    viable_path = [];
  end
end
